Computation of vertical fluid mobility of CO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}2, methane, hydrogen and hydrocarbons through sandstones and carbonates

Over the last decade, there has been an irreversible shift from hydrocarbon exploration towards carbon storage, low-carbon energy generation and hydrogen exploration. Whilst basin modelling techniques may be used to predict the migration of hydrocarbons through sedimentary basins on geological timescales, there remains little understanding of how fluids behave at the basin scale on present-day timescales. We apply the Darcy flow equation to present an algorithm to determine the basin-scale mobilities and maximum vertical velocity, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{max}$$\end{document}vmax, of CO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}2, methane, hydrogen and hydrocarbons with depth for sandstone and carbonate. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{max}$$\end{document}vmax for CO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}2 and methane are on scales of m/year, whilst values for hydrocarbon fluids are an order of magnitude smaller than for other fluids. Our results indicate that the fluid mobility of subsurface CO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}2 may be sensitive to surface and near-surface temperature variations. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{max}$$\end{document}vmax for hydrogen is approximately 2–10 times greater than hydrocarbon fluids, yielding important consequences for the future use of basin modelling software for determining hydrogen migration for exploration and storage.

whereby ρ w = density of water, ρ f = density of fluid and g = gravity (9.80665 m/s 2 ) 7 . The mobility of a fluid is therefore a function of both fluid properties and rock type of the carrier, whereas buoyancy is independent of rock type. The terms of Eq. (1) include several important quantities that must be evaluated separately using the appropriate relationships that describe fluid volume, viscosity, density and buoyancy. A detailed description of the Darcy flow equation is described in "Methods" section.
Two-and three-dimensional basin models may take several days to compute on a home computer. This could be improved by assigning a fluid velocity that is dependent on properties of the carrier rock and depth. However, determining the velocities of fluids through different rock types remains a significant challenge. The overwhelming majority of research on fluid flow through rocks is focused on laboratory studies of specific examples or case studies or unsuccessful application of basin modelling software to simulate real-time fluid flow 8 . As we shall see, the determination of fluid velocities in different rock types as function of depth will be of immense value to both basin and future modelling studies for hydrocarbon exploration, CCS and potentially green hydrogen exploration.
In this study, we apply various Equations of State and bulk properties of common rock types to different fluids encountered within the subsurface to calculate the variation of physical properties (e.g. volume and viscosity) with depth and surface temperature. Depositional parameter values for two generalised rock types, sandstone and carbonate, are applied in this study. For carbonate, parameter values for micrite are used. Depositional porosity ( φ 0 ) values are 0.41 and 0.51 and compaction wavelengths are 3.22 km and 1.92 km for sandstone and carbonate, respectively 7 .
We present an algorithm to calculate maximum vertical velocity, v max , as a function of fluid mobility and buoyancy using the Darcy flow equation (Eq. 1). The results of this study may have important implications for the application of basin modelling techniques to understand migration of non-hydrocarbon fluids and for threshold depths, below which fluids have low mobilities.
Using Athy's relation and the values in Tables 1 and 2, we use the following algorithm to calculate the vertical mobility and buoyancy of CO 2 and methane (outlined in Fig. 1): (1) Solve chemical equation of state (EoS) to calculate fluid volume from pressure and temperature (2) Calculate rock porosity (φ) and permeability at depth (3) Calculate fluid-water relative permeability (k rp ) and connate water saturation (S w ) from porosity (4) Calculate fluid viscosity from pressure and temperature (5) Calculate vertical fluid mobility and buoyancy Table 1. Component properties of hydrocarbon gases, hydrogen, CO 2 and water 7,[9][10][11] . M W = molecular weight, T C = critical temperature, P c = critical pressure, v c = critical volume, ω = acentric factor. Mw (g/mol) T c ( • K) P c (MPa) v c (m 3 /kmol) ω CO 2 44.010 304. 19  www.nature.com/scientificreports/

Results
The results for fluid mobility and buoyancy calculations for CO 2 and methane at normal geological conditions for surface temperatures of 0 • C and 20 • C are shown on Fig. 2. Gas phase mobility increases exponentially as depth decreases, with maximum mobility at the surface. Non-gas phase mobilities for both CO 2 and methane in both sandstones and carbonates are significantly lower than corresponding gas phase mobilities and decrease exponentially with depth. Fluid buoyancy follows a similar pattern. The choice of EoS has a minimal impact on calculated fluid mobility and buoyancy for CO 2 and methane. Fluid buoyancies increase exponentially upon the transition to gas phase. Methane mobility varies little with T s , with results indistinguishable for T s = 0 • C and 20 • C. Calculated vertical fluid mobilities for CO 2 are particularly sensitive to T s . Fluid mobilities for both CO 2 and methane are approximately double in sandstone compared to carbonate. Whilst the magnitude of vertical mobilities of fluids vary little with surface temperature, the depth of liquid-gas phase transitions vary considerably where critical properties are close to surface conditions (Fig. 7).
The equations used to relate relative permeability and normalised water saturation (Eq. 13) are designed for application to fluid systems with water, oil-water and gas-oil components. However, it is intriguing to extend this analysis to other fluids. The properties of several other important hydrocarbon fluids and hydrogen are listed in Table 1. We apply the algorithm proposed in this study to calculate vertical fluid mobility and buoyancy for hydrogen, dry/wet gas and various hydrocarbon fluids (Figs. 3 and 4).
Unlike for CO 2 , the vertical fluid mobilities and buoyancies of all other fluids considered in this study do not show significant sensitivity to T s . The magnitude of vertical fluid mobilities of hydrogen are up to twice those of dry/wet gas and over an order of magnitude greater than other hydrocarbon fluids (volatile, light and black oil). The buoyancy of hydrogen varies little with depth as all conditions within the Earth are beyond its critical parameters. Whilst above 0.36 km, the buoyancy of hydrogen is close to that of other gas phase fluids (methane, dry and wet gas, Fig. 4) and is > 8000 kg m −2 S −2 . Hence, this raises the important consideration that hydrogen velocities may be 2-10 times greater than those for other fluids and must be accounted for when evaluating hydrogen migration whilst applying conventional fluid flow modelling techniques suited for hydrocarbon migration. The vertical fluid mobilities for volatile, light and black oils are an order of magnitude smaller than those for other fluids, reflecting the fact that these remain in the liquid phase at all levels below the surface. Whilst results for dry gas are close to those for methane, results for wet gas resemble the pattern seen for CO 2 when T s = 20 • C. This is expected, as conditions within the Earth exceed the critical properties of wet gas (Table 1).

Discussion
The results of this study have several important implications. The velocity of fluids through rocks may be calculated directly by multiplying vertical fluid mobility and buoyancy, according to Eq. (1). Results indicate maximum vertical gas-phase velocities of 1.5-4 m/year (sandstone) and 1-2.5 m/year (carbonate) for CO 2 and methane, respectively. As both vertical fluid mobility and buoyancy decrease exponentially below the gas-supercritical or gas-liquid phase transition, our results indicate a minimum depth for CO 2 storage 0.36 km and 0.59 km when T s = 0 • C and 20 • C, respectively.
Basin modelling. Modelling of CO 2 and H 2 migration are dependent on understanding the present-day velocities of fluids. However, there exists no direct method for calculating fluid velocities as a function of lithology and depth. Attempts to apply basin modelling software to predict the effect of fluid injection on basin-scale pressure evolution have been largely unsuccessful 8 . Whilst the assumptions used in this work are relevant for water, oil or gas fluid systems in generalised sandstones and carbonates, our results could have important implications for future investigations for carbon storage and hydrogen exploration.
Sedimentary basins around the world are often characterised by complex geological histories. The PT-path of sedimentary basins will likely cross several regimes (e.g. different lines on Fig. 5). The methodology outlined in this study may be used to evaluate fluid velocities through different geological regimes (e.g. moving from normal conditions to regions of overpressure) by choosing the appropriate PT-path for discrete depth intervals. The methodology presented in this study, if combined with a directional vector, may provide a fast and computationally inexpensive means to calculate fluid velocities as a function of depth and angle of incidence between fluid particles and geological horizons, e.g. impermeable seals. A potential method to calculate the direction of fluid flow beneath geological seals is flowpath bending, which is widely used to model petroleum flow in basin modelling software. A planar seal with dipping angle α in the y-direction and lateral water flow with angle β to the x-axis and a hydraulic head with a dipping angle of γ are used to define the angle ψ , which indicates the direction of fluid flow as: whereby �ρ = difference in density between fluid and water and ρ w = water density 7 . When water flow is also assumed to follow the dipping of the seal, this becomes: Combining calculations for v max and directional vectors could significantly reduce fluid flow simulation processing time compared to leading basin modelling software packages that rely on the creation of nodal meshes to process fluid particle-rock interactions.  www.nature.com/scientificreports/ The range of uncertainties for S w becomes large beneath depths of 2 km (Fig. 9). Future work should focus on expanding our understanding about the relationship between porosity and water saturation in different rock types. New experimental data for a wide range of rock types will improve Eq. (14) and broaden the possibilities for the methodology described in this study.
It is widely accepted that hydrocarbons have the capability to migrate long distance through carriers on geological timescales. For example, Hantschel and Kauerauf 7 estimate the velocity of oil with density = 400 kg m −3 The results of this work indicate that CO 2 , methane, dry gas and wet gas may migrate several metres vertically per year ( > 1000 km/Ma) when in the gas phase. It is very encouraging that these values are consistent in magnitude with reports in the literature of hydrocarbon velocities of up to 1000 km/Ma, and indicate that the approximations used in this method are appropriate 20 . It is also encouraging that exponential increases in fluid velocities are concomitant with liquid-gas phase transitions given that this is independent of rock type. Furthermore, the magnitudes of vertical mobilities match those reported for hydrocarbons of fluids of similar composition and molecular weights to those analysed in this study (excluding H 2 ). Whilst the Darcy flow equation (Eq. 1) is widely accepted as the leading method to describe fluid flow in porous media, there are several important drawbacks. Averaging of microscopic features, such as discontinuities, holes, fractures and microporosities allows for upsaling to macroscopic length scales. However, a disadvantage of this method is the resulting low spatial resolution which constricts modelling of migration channels. Basin models typically cover large areas ( 10 × 10 km to 1000 × 1000 km) and are gridded and layered. Up to 500 grid points and 50 layers may be specified, depending on model complexity. Volume elements within each layer contain a constant facies bulk continuum approximation and treat all elements within a volume as constant 7 . Upscaling of physical properties from the core to grid size is necessary to capture the effects of physical features  Table 1 calculated using PR78, SRK and RK EoS for dry/wet gas and hydrogen and using appropriate EoS for hydrocarbon fluids [12][13][14][15][16][17][18][19] . Wet gas exhibits similar properties to CO 2 at T s = 20 • C, whilst dry gas exhibits similar properties to methane. Hydrocarbon fluid buoyancies are significantly lower than dry gas and methane and decrease with molecular weight. Hydrogen buoyancy varies little with depth. www.nature.com/scientificreports/ characteristic of specific lithologies. This is accounted for in part by applying appropriate upscaling factors. However, upscaling parameters may be amended as necessary for different rock types (e.g. using laboratory data) to improve overall model robustness in future studies.
Carbon storage and methane. The geological storage of CO 2 is dependent on the existence conditions that keep CO 2 in the liquid or supercirical phase and has been a subject of significant interest over the last three decades. Long-term storage of CO 2 under the ground at the subsurface is considered feasible for three types of geological formations: hydrocarbon reservoirs (depleted), deep saline aquifers and unminable coal formations. When considering potential sites for CO 2 storage, various geological characteristics must be taken into account, including characterising cap rocks to determine the effectiveness of a seal, if there are any abandoned or active wells which can compromise seal integrity and whether there is a sufficiently large and permeable storage formation 21,22 . With low levels of radiogenic heat production, shallow sandstone and carbonate layers within onshore sedimentary basins will typically exhibit lower pressures and temperatures compared to offshore depleted hydrocarbon reservoirs. However, the nature and the effects of surface and near-surface temperature variations on the physical properties of subsurface pore fluid-CO 2 systems remain relatively unexplored. CO 2 phase behaviour is primarily controlled by subsurface temperature and pressure conditions. Time series of subsurface temperatures and records of surface air temperatures indicate that subsurface temperatures are attenuated by several degrees up to depths of 10 m. However, intrinsic rocks are not influenced beyond this depth. For example, in warm climates the temperature 1 m below the ground surface remains below 30 • C whilst the ground surface temperature may reach 65 • C 23,24 . However, fluids transported via pathways created by heterogeneities intrinsic to rock volumes and fractures (e.g. faults), may play an important role in transporting heat from the surface to rocks at depth 25 . Low-frequency temperature signals (e.g. decadal climate change) may be retained at depth ( > 100 m) while high-frequency temperature signals (e.g. monthly) are generally retained at shallow depths ( < 10 m) 26 . Recent research indicates that fluids migrating along fracture networks through rock volumes may be significant propagators of heat, and the combined effects of urbanisation and global warming may reach more than 100 m below the surface 27-29 . Taniguchi et al. 28 measure the depth of deviation from regional geothermal gradients in several major Asian cities, with deepest occurring in Tokyo, at 140 m. Furthermore, recent investigations into the effects of large-magnitude urban heat islands indicate the creation of significant downward components of conductive heat flow in the shallow subsurface, which are supplemented by downward heat transport by groundwater movement 29 . Hence, elevated groundwater temperatures may have potential consequences on the state of CO 2 hundreds of metres beneath Earth's surface. As shown in Fig. 7, the pressure required to cross the CO 2 supercritical phase boundary increases with temperature. With a lack of corresponding increases in pressure, increases in temperature due to heat transportation from groundwater may cause CO 2 at depth to remain in the gaseous phase. In the absence of elevated subsurface temperature (e.g., from urban heat islands, igneous activity etc.), our results indicate that the depth of the CO 2 liquid-gas phase boundary changes from 0.36 to 0.59 km between surface temperatures of 0 • C and 20 • C, respectively. We interpret this to mean that CO 2 remains stable in a liquid form at shallower depths beneath regions of lower surface temperature. Therefore, it is our opinion that the long-term efficiency of CO 2 injection will be maximised at northern latitudes in regions not affected by urbanisation, igneous activity or other sources of increased subsurface heat flow. It is intriguing that changes surface and near-surface temperatures due to urbanisation and climate change over the last several centuries may affect the geological storage potential for CO 2 , and further research is required in this area.
The application of PT paths (Fig. 5) and EoS to determine fluid mobilities assume steady-state basin conditions and is appropriate for modelling past fluid migration during basin evolution, e.g., basin modelling. However, during CO 2 injection, forcing terms (i.e., inflow and outflow) must be considered to evaluate fluid properties at the injection site and surroundings. Whilst it is beyond the scope of this study to consider these in our algorithm, it is reasonable to imply that the algorithm proposed may be of use to model CO 2 migration within geological basins post-injection, over both human and geological timescales.
One limitation is that the cubic EoS applied in this study neglect fluid dissolution in water. It is commonly assumed that hydrocarbons do not dissolve in water, with the exception of methane. In the subsurface, methane dissolves in water only in the immediate neighbourhood of hydrocarbon pathways. Methane dissolution can be treated separately from hydrocarbons and is dependent on the amount of methane, pressure and temperature 30 . Another important consideration is that the physical properties of CO 2 differ from typical hydrocarbons. CO 2 is slightly polar and and dissolves well in water. Recent adaptations of cubic EoS (e.g., PR, SRK) have improved the accuracy of CO 2 -CH 4 -H 2 O systems significantly by including an association term named the Cubic Plus Association, or CPA factor. Results from mixing experiments of CO 2 , CH 4 and H 2 O gases, CO 2 -brine and CO 2 -acid gas-H 2 O systems show that in general, water content increases isobarically with temperature, decreases isothermally with pressure and increases with increased CO 2 concentration in the feed gas [31][32][33][34] . To explore the effects of surface and near-surface temperatures on groundwater-CO 2 and pore fluid-CO 2 mixtures in the subsurface, future research should focus on incorporating adaptations of cubic EoS and the application of rock properties for specific lithologies into the fluid mobility algorithm proposed in this study.
Gilmore et al. 35 evaluate the case of a fault cutting through three vertically stacked aquifers and seals based on data from a naturally CO 2 -charged aquifer at Green River, Utah. CO 2 has been escaping at the site along fault zones for several hundred thousand years. Their analyses indicates that despite leakage, the majority of trapped CO 2 remains stable when the fault permeability is is comparable or less than reservoir permeability. When fault permeability exceeds reservoir permeability, > 50% of injected CO 2 remains trapped after 1000 years. In reality, the ability for CO 2 and other fluids to migrate along faults depends on the relative differences in permeabilities between layers and faults, and will be determined by the geological complexity of the area. To determine migration pathways of CO 2 through multiple layers in geological basins, calculations of fluid mobility www.nature.com/scientificreports/ must be combined with adequate evaluation of carrier and fault permeabilities. Whilst the presence of impermeable geological seals and trapping mechanisms are of paramount importance during selection of CO 2 storage sites, risk of undetected fluid pathways capable of transporting groundwater into sedimentary rocks cannot be ignored. Examples of complex and vast systems of small-scale polygonal faults are likely to be responsible for fluid migration in post-rift sedimentary successions that were previously thought to be absent of migration pathways, e.g. Qiongdongnan Basin, China 36 . We believe that our algorithm offers an important starting point to evaluate CO 2 mobility (and potentially velocity) using a computationally inexpensive method which accounts for gas-liquid and gas-supercritical phase changes with depth. The principles of the algorithm proposed in this study also offers a potential means to bridge the knowledge gap between modelling of CO 2 plume migration, rock properties and including the crucial step of an EoS which may be modified to account for polar fluid-water interactions. It is our hope that this initial study leads on to further research in this area and development of methods capable of calculating the effects of surface temperatures on CO 2 velocities during and after CO 2 injection for storage purposes.
Hydrogen exploration. Momentum behind hydrogen exploration is growing, with jurisdictions such as South Australia granting or receiving applications for 18 exploration licenses by six different companies searching for natural hydrogen since February 2021 37 . Recent studies link the geological stratigraphic accumulation of hydrogen to the presence of shallow ( < 10 km) mantle rocks multi overlaid doleritic sills and aquifers 38,39 . Lefeuvre et al. 38 show that major faults may act as conduits for natural hydrogen and provide a migration pathway from shallow mantle rocks under the Mauléon Basin, Western Pyrenees, to the surface. An example of a major natural hydrogen discovery and implementation of reservoirs is in Mali, west Africa. In 2012, natural hydrogen was discovered and eventually connected to a fuel cell to supply electricity to the town of Bourabougu 37,39 .
Hydrogen is generally considered insoluble in water below pressures of 40 MPa, with the Peng-Robinson and SRK EoS decreasing in accuracy as pressure increases 40 . When the hydrostatic pressure gradient for any region is approximately ∼ 10.5 kPa/m, it is known as the normal pressure gradient. Abnormally high hydrostatic pressure gradients of ∼ 21.5 kPa/m have been encountered, e.g., in geopressured/geothermal zones, Gulf of Mexico, Niger Delta and North Sea 41,42 . Abnormally low hydrostatic pressure gradients < 10.5 kPa/m have also been encountered (e.g. Pennsylvania, Oklahoma 40 ). Using values for the normal pressure gradient, it is reasonable to assume that the EoS applied in this study are accurate for modelling insoluble hydrogen properties to a depth of ∼ 1760 m in most circumstances.
Clearly, the application of well known basin modelling techniques presents an attractive opportunity for hydrogen exploration. Application of the algorithm proposed in this study to model hydrogen migration at greater depths or in regions of high pressure may be possible with future research that focuses on amending widely used EoS to account for hydrogen solubility at pressures > 40 MPa. Our results indicate that future application of basin modelling techniques and software for assessment of short-term basin-scale hydrogen storage should account for fluid mobilities that are 2-10 times greater than those of hydrocarbons at the same geological conditions above depths of ∼ 1760 m where EoS remain stable.

Conclusions
We present an algorithm to calculate vertical fluid mobility and buoyancy for CO 2 and methane as a function of depth for generic sandstone and carbonate by combining solutions to chemical equations of state, viscosity and permeability calculations and the Darcy-Flow equation. Vertical fluid mobility and buoyancy in both sandstone and carbonate decrease exponentially with depth and are significantly greater for gas phases compared to liquid, dense liquid or supercritical phases. The depth of the phase transitions for CO 2 is sensitive to surface temperature, whereby gas-liquid and gas-supercritical transitions occur at 0.36 km for T s = 0 • C and 0.59 km for T s = 20 • C, respectively. Surface temperatures > 0 • C push the pressure-temperature profile of CO 2 into the gas-supercritical region. Hydrogen velocities may be approximately 2-10 times greater than those of CO 2 , methane, dry and wet gas and over an order of magnitude greater than for other hydrocarbon fluids. Dry and wet gas follow similar trajectories as methane and CO 2 . Vertical fluid mobility and buoyancy of other hydrocarbon fluids (volatile, light and black oil) are an order of magnitude smaller than those of CO 2 , methane and dry/wet gas.

Methods
Darcy flow equation. As described in Hantschel and Kauerauf 7 , the driving forces for fluid flow are pressure potential differences, where the fluid pressure potential is the pressure reduced by the pressure of a static fluid column with a corresponding vertical pressure gradient. The potential, u p for any phase p is thus defined as: with ρ p = density of fluid phase p, g = gravitational acceleration and z = depth. Darcy's law states that a potential difference causes flow according to: Here, v p = velocity of phase p and µ p = mobility. The mobility of in multi-phase fluid systems is usually split into three factors, such that: www.nature.com/scientificreports/ where k rock = rock permeability, k rp = relative permeability and ν p = viscosity of phase p. This provides a basic formulation of Darcy's law, with a comprehensive formulation given by: The driving force −∇u p is a gradient, which points in the direction of the steepest decrease of the potential field u p . It is multiplied with a tensor µ p ∝ k which describes the anisotropy of the rock permeability so that the resulting flow velocity v p is not necessarily pointing in the same direction as −∇u p . The Darcy flow equation may therefore be evaluated as: In the case of fluid flow in the subsurface, the driving force |∇u| is buoyancy between fluid phase p and water that is given by: Fluid properties in rocks. Fluids may be grouped into compounds with approximately the same physical properties, and are typically 'pure' (e.g,. CO 2 , methane or water) or 'lumped' chemical species (e.g., alkanes) 9 .
The subject of fluid analysis can be subdivided into three parts, the determination of the coexisting phases, their compositions and their properties. In this study, we consider both pure fluids (CO 2 , methane, and hydrogen) and lumped fluids (dry gas, wet gas and hydrocarbon liquids). The term 'phase' refers to the chemical state of a fluid (e.g., gas, liquid or supercritical). In large-scale modelling of fluid migration, it is usually assumed that a water phase is present. This is true for basin modelling techniques, where the water phase is commonly separated from hydrocarbon phases due to the non-polar nature of hydrocarbon fluid molecules. We extend application of these EoS to CO 2 , methane and hydrogen as the PR78 and SRK EoS only become unstable when the reduced temperature, T r < 0.6 , whereby T r = T/T c , T = temperature and T c = critical temperature 44 . As T c values of CO 2 , methane and hydrogen are 305.1 K, 190.6 K and 33.2 K, respectively, the corresponding minimum T r values for the PR78 and SRK EoS to remain stable are 183 K for CO 2 , 114 K for methane and 20 K for hydrogen. Hence, the PR78 and SRK EoS remain accurate within the range 0 • C < T s < 20 • C and at temperatures exhibited in most geological basins across the world. Other EoS that may be used include the Reidlich-Kwong (RK) EoS, which is suited for small, non-polar and short-chained molecules 12,13,15,19 . Other EoS have been developed specifically for use with hydrocarbons, e.g. volatile, light and black oils 17,18 . Component properties for CO 2 , methane, hydrogen, dry gas, wet gas and hydrocarbon fluids are listed on Table 1. Whilst EoS may be applied to calculate fluid pressure, temperature and volume, viscosity is another important indicator for phase property characterisations and fluid velocities within rocks and geological basins. Viscosity is often modelled as a quantity dependent only on pressure, temperature, density and the amount of dissolved gas 9 . Whilst advanced theories, such as the friction-theory or free-volume model match laboratory data well, a lack of field data and unknown component parameters (especially for heavy hydrocarbon compounds) make these unsuitable for basin modelling 7 .
Compared to methods that require laboratory or field data, direct approaches such as the empirical Lohrenz-Bray-Clark (LBC) model are simpler to apply in basin modelling scenarios where data is often scarce 12,13,45 . Viscosities can be evaluated very fast due to the simple nature of the LBC-formulas. However, models with lower performance are often not usable in fluid flow simulators. It must also be noted that the LBCmodel is based on a polynomial of degree 16. Polynomials of such a high degree are known to easily become numerically unstable and therefore LBC-based models must be evaluated with care 7 . Furthermore, the pressure correction for gases in the LBC method had a tabular formula and was not presented entirely 45 . Using this distinction introduces a discontinuity between the liquid and gas viscosity. A modification to the LBC method may be made by applying the Herning mixing rule 46 . This was validated to successfully model multicomponent and multiphase fluid viscosity using the modified LBC approach for an example using the independent Stiel-Thodos method 13,47 . Hence, this modified LBC method is preferred over the original LBC method 13,45 . www.nature.com/scientificreports/ The range of pressure and temperature (PT) paths in in Fig. 5 defines an area of possible PT points in arbitrary sedimentary basins. To simulate generic geological conditions, a geothermal gradient of 25 • C/km and PT gradient for normal conditions (0.5 MPa/K) was applied (Fig. 5) to calculate the molar volume and density of CO 2 and methane by solving the PR78, SRK and RK EoS 7,12,13 . The modified-LBC method was used to calculate fluid viscosity (Fig. 6 13 ).
The transition between phases (solid, dense liquid, liquid, gas and supercritical fluid) is controlled by a fluid's pressure and temperature relative to its critical pressure, P c , and temperature, T c . Figure 7 shows phase diagrams of CO 2 and methane with PT paths for T s = 0 • C and 20 • C. For CO 2 , P c = 7.3773 MPa and T c = 30.98 • C and is encountered at 0.59 km under normal geological conditions. When T s < 20 • C, CO 2 undergoes three phase transitions at 0.36 km (gas-liquid), 0.59 km (liquid-dense liquid) and 1.24 km (dense liquid-supercritical). However, when T s ≥ 20 • C, only one phase transition occurs at 0.59 km (gas-supercritical). For methane, all  [14][15][16]19 . Blue dashed lines represent supercritical-dense liquid and dense liquid-liquid phase transitions for T s = 0 • C. Black dashed lines represent gas-supercritical phase transition for T s = 20 • C (see Fig. 7). The change from liquid to gas for all fluids is marked by an exponential increase in molar volume and decreases in viscosity and density.  (Table 1) and rocks (Fig. 8). Another important parameter is porosity, φ , and is defined as the ratio of free pore space to rock volume. φ controls the volume between grains in a sedimentary rock which may hold fluids and reduces as with compaction. The variation of φ with depth must be considered when evaluating fluid velocities through rocks. Athy's relation parametrises φ using: whereby φ 0 is the maximum (depositional) porosity and k is the Athy compaction parameter (Table 2 48 ). The flow of fluids through rocks generally involves more than a single phase. The ability of one fluid to flow through a rock or reservoir is affected by the presence of other fluids. Relative permeability, k rp describes multiphase flow in reservoirs as the ratio of the effective permeability of a fluid to the absolute permeability of the rock, k rock . The effective permeability is a relative measure of conductance of the porous medium for one fluid phase in the presence of other fluid phases.
Approximating k rock by applying various relationships between porosity, φ and log(k rock ) is well established and may be performed using several methods and applying appropriate upscaling. These include the multipoint method 7 or revised Kozeny-Carman relation 49 for practical use in basin modelling. The multipoint method describes the relationship between φ and log(k rock ) as a piecewise linear function (Fig. 8). A compilation of porosity-permeability data points for different lithologies are tabulated in Appendix A of Hantschel and Kauerauf 7 . The revised Kozeny-Carman relation 49 describes k rock as a function of porosity, lithology-dependent specific area and scaling factor. The multipoint and revised Kozeny-Carman relation yield similar values for k rock (φ) (see Figure 2.16 in Hantschel and Kauerauf 7 ), however as the latter is not suitable for carbonate rocks we restrict our analyses to the multipoint method. Combination with Eq. (12) allows computation of k rock as function of depth and lithology-dependent initial porosity and compaction parameters.
Vertical and horizontal rock permeabilities for use in basin modelling may be calculated from hand-specimen measurements by using appropriate anisotropy and upscaling factors. Anisotropy is the ratio of horizontal and vertical permeability and is dependent on rock type. Basin scale values for horizontal and vertical permeabilities are calculated from hand specimen values multiplied a horizontal and vertical upscaling factor, respectively. The resulting higher values for larger scales are caused by macro-fractures, inhomogeneities and permeable inclusions 7 . Here, we apply typical upscaling factors of 1 (vertical) and 50 (horizontal). However, greater upscaling factors (e.g. vertical = 10 and horizontal = 500) to basin scale elements with lengths greater > 50 m are reported in the literature for sandstones 50 .
Numerous models for estimating relative permeability have been proposed over the last half century. Most relative permeability models assume relative permeability to be a function of connate water saturation ( S w ) and connate gas saturation ( S g ) respectively 51 . The relative permeability of any fluid is zero below its critical saturation, where it becomes immobile. Saturations are for often rescaled into normalised or effective saturations ( S e ) which map the saturation interval between the connate and the critical saturation to an interval of 0 < S e < 1 . The relationships between normalised saturations of water, S we , gas, S ge and oil-gas, S goe , are as follows: whereby k rw , k row , k rog and k rg are the relative permeabilities of water, oil-water, oil-gas and gas components respectively. k rw and k rog represent water-liquid flow whilst k rg and k rog represent liquid-vapour flow, respectively.
for k rw , k row S goe = S g /(1 − S wc ) for k rog S ge = (S g − S gc )/(1 − S wc − S gc ) for k rg  www.nature.com/scientificreports/ and k rog = f (S g ) . Assuming that phases do not interact during flow, the flow of each phase can be treated as if the other phases are part of the solid rock matrix 51 . In this case, the relative permeability of oil, k ro = k row k rog . The quadratic relationship of Hantschel and Kauerauf 7 may therefore be applied to determine relative permeabilities from normalised effective saturation values. We use this relationship to approximate the flow of the fluids (Table 1) in sandstones and carbonates.
Water saturation from porosity. Determining water saturation, S w , is an extremely challenging petrophysical calculation to perform. S w is used to quantify the hydrocarbon saturation, (1 − S w ) and must be evaluated in basin modelling. Complexities arise because there are a number of independent approaches that can be used to calculate S w . In wellbores, the following methods can be used to determine S w : • S w calulations from resistivity logs and by application of a model relating S w to porosity, φ , connate water resistivity and lithology-dependent electrical properties. • S w calculations from laboratory capillary pressure and saturation ( P cap /S w ) measurements by application of a model relating S w to various rock, fluid properties and height above the free-water level. • S w calculations using oil-based mud (OBM)-core-plug Dean-Stark-water-volume determinations.
• Combinations of these methods The amount of data that is available often dictates which method to use to determine S w . For the purposes of basin modelling, which requires the application of generalised lithologies and fluid types, S w must be calculated from data input by the user, e.g. rock type, fluid type and depth. Laboratory-based methods for measuring S w are not suitable for this work, as we are concerned with the properties of generalised rock types and fluids on the basin-scale.
The results of early experiments indicated that the product of φ and S w is constant 52 . The magnitude of this constant was shown to be a related to rock type and indirectly to permeability, k. Better quality rocks were found to correspond with low constant values. Extensive analysis of core data and petrophysical estimates of porosity and irreducible water saturation, from all types of reservoirs worldwide, suggests that this relation is a unique solution to a more general equation: whereby the value of the power function Q ranges from 0.8 to 1.3, with many reservoirs close to 1.0. For sandstones, 0.02 < C < 0.10 and 0.005 < C < 0.06 for carbonates. The values of Q and constant C can be easily derived by plotting logφ against logS w . The gradient of this linear plot = Q . Projection of the straight line against φ = 1.0 gives the value of the constant, C 52,53 . Figure 9 shows the relationship between S w and φ calculated using Eq. (14) and depth from Eq. (12). Uncertainties may be calculated using the upper and lower bounds of experimental constants and exponent terms from laboratory measurement 53 . Uncertainty increases significantly as porosity is reduced through compaction ( φ < 5% ) and with depth. Uncertainties increase significantly below 2 km and cover almost the entire range of possible S w values by 4 km depth. Hence, we recommend that application of this methodology should be limited to depths above 2 km. Future laboratory measurements of coefficients of Eq. (14) may lower uncertainties and increase model robustness beneath 2 km. (14) φ Q × S w = C Figure 9. Connate water saturation ( S w ) as a function of porosity and depth calculated using parameters for sandstone and carbonate and Athy's relation 7,48,53 . (A,B )= sandstone, (C,D) = carbonate. Grey regions represent uncertainties calculated using upper and lower bounds of experimental constants for Eq. (14) 53 . Uncertainties beyond depths of 2 km become significant, hence we do not extend this method beyond this depth.